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Abstract 

We want to compute the cumulative distribution function of a one-dimensional Poisson stochastic 

r 

integral I{g) = / g(s)N{ds), where A'^ is a Poisson random measure with control measure n and g is 
Jo 

a suitable kernel function. We do so by combining a Kolmogorov-Feller equation with a finite-difference 
scheme. We provide the rate of convergence of our numerical scheme and illustrate our method on a 
number of examples. The software used to implement the procedure is available on demand and we 
demonstrate its use in the paper. 



1 Introduction 

Let T > 0, and let N{-) be a Poisson random measure defined on the interval [0,T], with the Borel cr-field, 
and control measure n{ds), which we assume to have a density n{ds) — n{s)ds. For appropriate functions 
g{s), s G [0,T]. the Poisson stochastic integral 

1(9) = I 9{s)N{ds), (1) 



is a random variable defined as the limit in probability: 

71 

P- lim y"Q{s^)N{[s„s,+i)), A„ = {0 < si < S2 < ••• < s„ < T}. (2) 
For this limit to exist, we require that kernel g satisfies 

min((7(s), l)n{ds) < oo. (3) 
The characteristic function of /(g) is expressible in terms of the control measure and kernel, and is given by 

(j)gie) = Ee'^-f^s) ^ exp ( / (e^'^s^"' - l)n(ds) ) . (4) 
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For more details about Poisson integrals and some discussion of applications, see Chapter 3. 

It is important to note that in the case where n is a finite measure on (0,T), the distribution of I{g) is 
that of a compoimd Poisson distribution. This can be seen in the case where g is strictly increasing and has 
inverse g^^ by making a change of variable u = g(s) in the integral in 

exp (^J^ (e*^»(") - l)n(ds)j = cxp (^J^ ^ \e*^" - l)n{u)dv^ 

where fi{u) = {g~^y{u)n{g~^{u)), which is the characteristic function of a compound Poisson distribution 
(see Proposition 1.2.11 in [5]). Thus, I{g) in this case has a compound Poisson distribution with rate 

r = jg-i(^Q-^ n{u)du — n{s)ds and jump distribution with density r~^n{u), g~'^{0) < u < g~^{T). A 
similar argument shows this for non-monotone g, except one must construct n carefully by breaking the 
integral ([T]) into pieces on which g is either flat or strictly monotone. Thus, by considering stochastic 
integrals of this type, we are covering a large class of interesting distributions. 

The distribution of integrals such as ([TJ would be easy to obtain if the measure was Gaussian. In this 
case I{g) is also Gaussian: I{g) ^ A(0,/g g'^{s)n{s)ds). Since N is Poisson, however, the integral I{g) is 
not in general Poisson and its distribution is not easy to determine. Our goal here is to study the cumulative 
distribution function (CDF) of I{g), which we denote by F{x) = Prob(/((?) < x) and develop a convenient 
numerical scheme to evaluate it. In particular, we focus on the following: 

1. A Kolmogorov- Feller type evolution equation associated to F. 

2. Smoothness properties of F. More specifically, we show that under certain assumptions on g and n, F 
lives in the Holder space C'^''^ {{g{Q), g{T)), for some < 7 < 1 depending on g. 

3. A numerical method for computing F. 

At first glance, one might attempt to compute the CDF F using the characteristic function ^ and the 
integration-based inversion methods of Abate and Whitt [Tj. For continuous F, their approach boils down 
to computing the following integral numerically: 

F{x) = - r Rei(l}g)iuf^^^du. (5) 

However, this method is not always efficient in this setting for multiple reasons. First, evaluating the 
characteristic function (/)g is not always easy, since /□'"(e*"^'-*-' — l)n{ds) is unlikely to have a closed form, 
and numerically computing this integral might be difficult. This adds another source of error on top of the 
"truncation" and "discretization" errors associated with the integration of ([5]). Moreover, these errors are 
difficult to bound exactly. Secondly, we will see in the following that F is not differentiable in general, which 
translates to a slow rate of decay of (j)g. This makes numerical integration difficult. 

We propose an alternative method for computing F which does not involve (directly) the characteristic 
function (|4|). Our method has several advantages over brute- force integration of ([5]). First, our method 
only requires knowing g and the density of the control measure n, and thus is much faster than integration 
methods which must evaluate the characteristic function (|4]). Second, we provide here error bounds given 
general assumptions on g and n. Third, our method generates the CDF F on an entire interval instead of 
a single point. Fourth, we obtain the more general cumulative distribution functions F(x,t), t > 0, of the 
following stochastic process: 

X{t) = Ao + / g{s)N{ds), 0<t<T, (6) 
Jo 
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where Xq is an independent random variable with given CDF Fq{x). In the sequel, we illustrate some of 
these advantages on an example. 

Our method involves solving an evolution equation satisfied by the CDF of X(t) with initial condition 
_Fo. A nice way to view such an evolution equation is to consider the characteristic function (j)g{uj,t) of ([5]) 
(we will from hereon drop the "g" subscript from (j)). Using we have at time < t < T, 

^{e,t) = 0o(0)exp (^V^'"*'^ - l)n{ds)^ , (7) 

where 0o is the characteristic function of Xq. Differentiating with respect to t gives 

^0(0,t) = ((e"^f(*)-l)n(t))0(0,t). (8) 

Thus, the characteristic function (j){uj,t) can be viewed as the solution at time t of the simple ordinary 
differential equation ([8]) with initial condition (f>{0,O) = 4'o{0). With this ODE perspective, the naive 
approach of brute-force integration can be thought of as a two step process: (i) evolving 0o under the 
trivial dynamics of ([S]) to obtain (p{d,t), and then (ii) computing ^ to get F. 

An alternative approach of obtaining an evolution equation for F is to notice that the process X in ([6]) is a 
continuous time Markov process, hence it satisfies the Kolmogorov- Feller forward equation ([3], Section 5.1). 
In the infinitesimal time interval [t,t + dt), X may take a jump of size g{t) with probability n[t)dt + o{dt), 
thus, the Kolmogorov-Feller equation for the single-time density p{x,t) of the process X takes the form 

t) ^ -n{t){p{x, t) - p{x - g{t), t)), t > 0, a; e M. (9) 

Integrating this equation suggests that a similar equation should hold for F, however some technicalities 
arise because p can contain "atoms". We show directly that a similar equation holds for the CDF, with 
the exception of points where the CDF is discontinuous. We then solve this evolution equation numerically 
using a finite-difference scheme. 

This paper is outlined as follows. In Section [2J we state a theorem which gives an evolution equation 
of the function F. We then study the smoothness of the function F{-,t) in Section [3] We then describe 
in Section |4] a numerical method for calculating F using a finite-difference scheme, and provide a rate of 
convergence for our method. In Section [SJ we apply our methods to a collection of examples and in Section 
[SI we establish lemmas related to the quality of the approximation. Finally, in Section [7] we provide a guide 
to the software, written in MATLAB, which allows the user to obtain the cumulative distribution function 
and density function of the integral g{s)N{ds) for general g. 



2 Kolmogorov-Feller equation for F 

In this section we show that the CDFs F{x,t), i > of X{t) satisfy a Kolmogorov- Feller-type equation. As 
mentioned in the introduction, the form of this equation is easy to guess based on the Kolmogorov-Feller 
equation © satisfied by the density of the process X{t). 

Since we will be dealing with distribution functions, we must be mindful of discontinuities which may 
arise. Given a CDF F, we will let C{F) denote the points of continuity of F. Also, we define the right-time 
derivative as 

^Fix,t)= hm Fi.,t + h) Fix,t)^ 

dt+ ^ ' h^o+ h ^ ' 

The following theorem specifies the equation for F which will be used the sequel. 
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Theorem 2.1 Let g{s), s > be a right- continuous function satisfying (0), and consider the process 

X{t) ^Xo+ f g{s)Nids), 0<t<T, (11) 



where N is a Poisson random measure with control measure n(ds) = n{s)ds and Xq is a given random 

variable. Let F[x,t) be the CDF of X at time t. Then, for all {x,t) such that x — g(t) € C'{F{-,t)), F 
satisfies the following equation: 

-^F{x, t) = -n{t){F{x, t) - F{x - g{t), t)). (12) 

Proof. Let x,t > such that x — g{t) e C{F) and let /i > 0. Consider the difference 

F{x, t + h)- F{x, t) = E[l(_^,,](X(t + h)) - l(_^,,](X(i))]. (13) 

By conditioning on the number TV of Poisson arrivals in the interval [t, t + h), the expectation on the right- 
hand side becomes 

oo 

E[l(_oo,.](^(i + h)) - l(_^,,](X(i))] = J2 E[l(-oo,..](^(t + h)) - (X(t))|iV ^ n] PToh[N = n] 



n=0 



E[l(-oo,.)(^(t + M)-l(-oo,.](^(i))|A^ = l] / n(r)rfr e-/.' "M''- + ©(/j^)^ (^4) 



t+h 



Here we have used that facts that ii N ^ 0, X{t) = X{t + h), and the probability of seeing two or more 
Poisson arrivals in the interval [t, t + h) is 0{h^) since we've assumed n is continuous. In the event of one 
arrival in [t, t + h), namely TV = 1, we have X{t + h) = X{t) + g{S), where S £ [t,t + h) is a. random time 
with density 

-fTTT^ , se[t,t + h) 

J^^'nir)dr . (15) 

0, else 

Conditioning on S, now becomes 

t+h \ I r-t+h \ 

]E[l(_oo,.](^(<)+.9(-'*))-l(-oo,.](^(i)) \ S = s\fs{s)ds \ n(r)drj "(^^'^^ + 0(/i2) 

t+h \ 

E[l(_^.,] {X(t) + gis)) - l(_oo..](X(0)]n(s)ds J e" "t^)"^- + 0{h^) 
IE[l(-oo,.-g(s)] (Xit)) - (X(i))]n(s)ds j e- "(^^''^ + 0(/i^) 

t+h \ 

{F{x - g{s), t) ~ F{x, t))n{s)ds + 0{h^). 
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Now, since F{-,t) is continuous at a; — g(t), g is right continuous, and n is continuous, let ft, — > 0^ and got 
Fix,t + h)-Fix,t) ^ ^.^^^ (h-' r\Fix-gis),t)~F{x,t))n{s)ds]e-f^'^""^-^''^ (16) 



= 7i(i)(F(a;-g(t),0-i^(a:,<)), (17) 

which completes the proof. 

■ 

Remarks. 



1. If F and g are continuous, then the right and left time derivatives of F coincide and the left hand side 



of (fT2)) can be replaced with MrFiXjt). 



2. If a; — g{t) ^ C{F) and 5 > is continuous and monotone increasing, then x — g{s) x — g{t) as s \t 
and the limit in ([T6| still exists, except we have instead 

^=n(0(F-(.T-.g(i),t)-n^,i)), (18) 

where F^ is the left-continuous version of F, i.e. F~{x,t) = limj^_j.2;- F{y,t). This fact will be used in 
the sequel. 

3. The continuity set C{F) can be identified from g and n as follows. We can write 

N 



I{g)=Xo+ f g{s)Nids)^Xo+Y,giS,), (19) 
Jo „_i 



where is a Poisson random variable with mean n{s)ds. Then, the atoms of I{g) correspond to 
the atoms of the random sums Xq + J^^^i gi^i)^ with n = 0, 1, 2, . . . and Si are the "Poisson events", 
which are i.i.d. with density 

0<s<t. (20) 

Jo n[T)dT 

The continuity set C{F) is the complement of this set of atoms. 

For example, if Xq = 0, _g = is a constant function, and n(s) = 1, then the atoms of the random 
sums X]i=i g{Si)j which is the set {nrj, n = 0, 1, 2, . . . }. 

4. To arrive at Theorem 12.11 we used properties of the Poisson integral to derive the evolution equation 
(|12[) . This reasoning can also be turned on its head - this theorem can give a probabilistic interpretation 
to initial value problems for a class of differential-difference equations of the form 

du 

— = -n{s)(u{x,t)-u{x-g{t),t)), xeM.,t>Q ^^i) 
u{x, 0) = ua{x) 

where n, g and are given functions and u is an unknown function of x and t. By computing the 
CDF of the Poisson stochastic integral ([1]), we are essentially computing the fundamental solution of 
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this differential-difference equation, since the solution to (|2ip can be written as the convolution of uq 
with the CDF F{x,t) of the integral ([J) up to time t: 



u{x,t)^ / uoix-u)dF{u,t). (22) 

— oo 

To illustrate this point, consider the following example: Let b{t), i > be a given smooth function 
with derivative b'{t), and recall the classical transport PDE: 

du , I , ,du ^ 

^ = -^(%' ^^^'^^0 , (23) 
u{x, 0) ~ uq{x) 

where, for simplicity, the initial condition uq is some smooth and bounded function. The solution to 
this equation is given by 

u[x,t)=uo[x-b{t)), t>0. (24) 

Now, let /i > 0, and say we want to interpret the solution to the following variation of ()23p where we 
replace du/dx with a finite-difference approximation: 

dt h , xfcK,r^u _ ^25) 

u{x, 0) = uq{x) 

Notice that (P5|) is of the form with g{t) = h, and n{t) = b'{t)/h. This corresponds to the trivial 
Poisson random integral 

' hN{ds), n{s) = (26) 
h 

which has the distribution of the random variable hN{[Q, t]), where A^([0, t]) has a Poisson distribution 
with mean b{t)h^^. 

Observe that in this case, (l22t can be expressed as u{x, t) — EuQ{x—hN{[0, t])) which is the convolution 
of uq with the CDF of hN{[0, 1]). Compare this to the solution of the classical transport equation p4|. 



3 Smoothness of F{-,t) 

Our goal is to develop a numerical scheme for approximating the CDF F which is based on the differential- 
difference equation (jl2p . In order to estimate the error in our approximation to the true CDF, we must 
have some smoothness properties for F, and thus make some assumptions on g. We will assume that g has 
a continuous inverse. While this seems like a restrictive assumption, we discuss ways of generalizing our 
method to non-invertible g in section (|4.3p . 

In the following, we assume that g is a non-negative, strictly monotone continuous function on [0, t] with 
inverse g~^. We will also suppose that n{s) is a bounded function. Since g > 0, the integral J* g{s)N{ds) 
is zero if there are no Poisson events in [0,t), i.e. if N[0,t) = which happens with probability e"^ 
Also, J^g{s)N{ds) > 5(0) if 7V[0,t) > 1. Therefore, F{x,t) = for x < 0, F has a jump of size e' J'>^'^'^' 
at x = 0, and F = e~-/o on the interval [0,g(0)) (see Figure [1]). Having identified the discontinuity 

at X = in the CDF of F, we now focus on the following question: How smooth is F for x > 0? This will 
depend on g. 
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ST 



Jl* n(s)ds 



{ 





X 



g(0) 



Figure 1: This shows the behavior of F{x, t) near a; = with t fixed. We see a jump of size exp(— n{s)ds) 
at a; = whieh equals the probabihty of zero Poisson events in time {0,t). Also, F{x,t) is eonstant in the 
interval [0,(7(0)) since with one or more Poisson event in (0,t), g{0) is the minimum value the integral can 
take. 

Let [/ C M. For functions u defined on U and for < 7 < 1, recall the definition of the Holder seminorm 
[•]co-7((7) and the Holder space C°''^{U) defined as 



Mco.T(C/) 



sup 



\u{x) - u{y)\ 



yeu \x-y\'' 



C^'^U) = {u : Mco..(c/) < 

The exponent 7 is called the Holder exponent of u. For more on these spaces, see [3] section 5.1. 
The next theorem shows that F{x, t) for a; > and t > Q lies in the same Holder space as 5"^. 



(27) 
(28) 



Theorem 3.1 Let g he a non-negative, strictly increasing continuous function on [0,7] with inverse g ^. 
For < t <T fixed, let F{x,t) — Prob[X(i) < x], where X{t) is defined by 



(29) 



X{t) = / g{s)N{ds), 

JQ 

where N is a Poisson random measure with control measure n{s)ds with n a bounded function. Let 



1* = sup n{s). 

0<s<t 



(30) 
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Then, ifg-^ e C°^'^{{g{Q),g{t))) with < 7 < 1, then F{-,t) G C"'^'^((0, 00)) and 
where the constant Cn,t is given by 



C„ 



(1 - e^o 

t 

lo nis)ds 



Proof. Fix t > and let S, Si, S2 ■ ■ ■ be i.i.d. random variables with density 

, , n(s) 
nt{s) = f , , . , 0<s<t. 



Jo n{T)d7 



Let nt{x) = nt{s)ds be the CDF of S and let 



Gix) = Prob[g(S') < a;] 



(31) 



(32) 



(33) 



'0 x<giO) 

Prob[S' < g-'^{x)] g{0) <x< g{t) 
1 x>g{t) 

'O x<g{Q) 

nt{g-\x)) g{0)<x<g{t) . 
1 x>g{t) 

By conditioning on the number of arrivals TV in the interval [0, t], we have for x,y > 0, 

00 / n n \ 

\F{x,t) - F{y,t)\ = J2 P™b[^g(50 < x] - Prob[^g(5.) < y] P„ 

n=l \ 4=1 1=1 / 

where P„ = Prob[A^ = n]. Let G*" denote the n-fold convolution of G with itself. Now we have 

00 

\F{x,t)-F{y,t)\ < \G{x)-G{y)\Pi + Y,\G*"{x)-G*"{y)\Pn 

= \G{x)-G{y)\Pi + y] / G{x - u)G*^''~^\du) - G{y - u)G*^"-^\du) 



(34) 



(35) 



< \G{x)-G{y)\Pi + J2 



71^2 

nmax{x,y) 



n=2 



\Gix -u)- G{y - u)|G*("-i)(dM)P„ 



Pn 



(36) 



Thus, for < 7 < 1, 



\F{x,t)~F{y,t)^ ^ \G{x)^G{y)\ 



,1=2 -^0 



max(a;,y) 



n=l 

= (l-e-/o"(^)^^)[G]co,.((o,oo)) 



M:^^if)_a^G^(-i)(d,)P„ (37) 

(38) 
(39) 



What remains is to show that [G]co,y((o^oo)) is bounded by a constant muhiple of [g ^]c'''^(s(o),g(t))- 
Observe that fi't{x) = nt{x) where nt is defined in (p3)) and thus |fi((a;)| < ■ {J^ n{s)ds)~^ , where = 
supq<;^<j n{s). In view of p4l) . the mean value theorem imphes that for g{0) < x,y < g{t), 

\G{x)-G{y)\ _ \nt{g-\x))-Mg-\y))\ n\ \g-\x) - g"\y)\ n\ , 

(40) 

If < X < git) < y, then G{y) = 1 = = nt{g-\git))) and 

IGjx) - G{y)\ ^ IMg-H^)) - ntjg-Hgmi ^ r -11 

F^M^ -J^s^' ho^nm^m- m 

A similar argument shows this also holds for x < g{0) < y. Finally, if a;, y > g(t) or x,y < g{0), then 

]cO.^(g(0),s(t))- 



G{x) - G{y) ^ < * , [g-'' 



Thus, in all cases, 

\G{x)-Giy)\ 



^ rt * , [g ^]c0'"'(g(0),3(t))- (42) 



/on(s)(is 
Hence, (I2H1) and (gl]) imply 

[^(■,i)]cO'"'((0,oo))^ , [g~^]c"^y{g{0),git)) = ^nd [ff"^]c«--'(ff(0),g(t)) ■ (43) 

Jo «(s)as 

which finishes the proof. 

■ 

In the preceding Theorem, we assumed g is monotone. What happens if g is piecewise monotone, i.e. g{t) 
is monotone on the intervals [ti,ti^i)^ for i — 0,1,..., M? Then, the integral of g over [^Oi^m) can be 
"constructed" inductively by considering the sums 

h^h-i+ l\{s)N{ds), i = l,2,...,M (44) 



with Iq = 0. Since A'^ is a Poisson measure, the summands in (j44p are independent random variables. Hence, 
the full CDF of g{s)N{ds) can be seen as the successive convolution of the CDFs of the /^'s above with 
the CDFs of the integrals J*^ ^ g{s)N{ds). Wc know the smoothness of the CDF of each integral, but what 
about the smoothness of a convolution of two such CDFs? The following corollary shows that the Holder 
exponent associated to the convolution of any two of these CDFs will, at worst, be the lesser of the two 
Holder exponents associated to the individual CDFs. 

Corollary 3.2 Let g, n and X(t) and F be as in Theorem h3.1\) . and suppose Xq is a non-negative random 
variable independent of X{t) with CDF Fq G C*'''''"(0, cxj). Then, the CDF F of the sum Xq + X{t) lies in 
the Holder space C'''^(0, oo), where 7 = min(7,7o), and 

[F]c»r, (0^00) +max(l, [-Fo]co.™(o,oo)), 7o > 7 

[F]co,,(o,oo) < ^ [F]cOn{o,oa) + [^o] C«-^o (0,oo) , 70 = 7 (45) 

[Fo]c'>.ya{o,oo) + max(l, [-F]c«.t(o,oo)), 7o < 7 
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Proof. Let < y < X. Observe first that if 7 > 7, then 

t"^^^"'^- \^-yP - l-y|<i 

I \x-yP 

< max(l,[F]co,.(o,oo)). (46) 

Similarly, if 70 > 7, wc have 

[Fo]c«-^ < max(l, [i^o]c»'^o(o,oo))- (47) 
To compute the Holder norm to F, we write it as the convolution of F and Fq: 



Fix)^F{y)= / {F{x^u)^F{y~u))dFoiu). (48) 
Jo 

Since = on the negative real axis, the above simplifies to 
_ _ fv 

F{x)-F{y)^ {F{x~u)^ F{y~u))dFo{u)+ F{x ~ u)dFo{u). (49) 



Now, Theorem p.ip and (|46p imply that 



\j^{F{x^u)-F{y-u))dFM ^ F \F{x-v)^F\y~u)\ 



\x-yV Jo \x-yr 



< , (50) 

I max(l, [i^]co,7(o,oo)) 7>7 



And, dH]) implies 



!yF{x^u)dFo{u) ^ Foix)~Foiy) ^ | [FqIcc™ (0,00) , 7o = 7 ^^^^ 
~ ~ [max(l, [Fo]co,-,o(o,oo)) 7o > 7 

Thus, (gH), together with ([50]) and ([511) imply the result. ■ 



CoroUarv 13 . 21 points to a possible "worst case" result about the smoothness of a convolution of two CDFs 
with different smoothness properties on (0, 00) by saying this convolution on {a; > 0} will lie in the larger of 
two Holder spaces {C^'^{0, 00) and C'^ '^''{0, 00)). This worst case occurs if F and Fq are both discontinuous 
at X = 0, as indicated in the following corollary. 

Corollary 3.3 Under the assumptions of Corollarv \3.2[ suppose further that Fq has a discontinuity at x = 
and that C^'"'{0, 00) and C^''^"(0, 00) are the smallest Holder spaces to contain F and Fq, respectively. Then, 
C°'^(0,oo) is the smallest Holder space to contain F. 



Proof. Notice that if Fq and F are discontinuous at x = 0, the quotients ((50)) and ([51]) can also be bounded 
below as 

\x-yp \x-yp 
ly F{x ~ u)dFoiu) ^ ^ Fo{x)-Fo{y) 
\x-yp - \x-yp 
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And thus it follows from that if C'^^'' and C'^'"'" arc the smallest Holder spaces to contain F and Fq, 
respectively, then taking supremums over all x and y above will cause one of these two ratios to become 
infinite if 7 > min(7,7o). Thus, (|48p implies that C*^'^ is the smallest Holder space to hold on {x > 0}. 

■ 

Notice that for sums like (|44|) . the corresponding CDFs will typically contain a jump at a; = 0, thus the 
"worst case" is actually what occurs. 



4 Obtaining F{x,t) 



There are standard finite-difference schemes for solving linear PDEs (see for instance [5], Chapter 9). In this 
section, we describe such a scheme which can be used to solve ([T2|) iteratively on the interval < t <T. We 
also study its convergence properties. 

To begin, we need to define various components of a discrete approximation of F(x, t) which we will 
denote F{jS,ih) and define as follows. Fix an interval [0, L] and subdivide it with a step size of 6. That is, 
consider 



where 



As = {j6 : 0<j <L/6} 

= {a;o < xi < a;2 • • • < xm}, 

M^\As\ =L/S+1 



(52) 



(53) 



is the size of the mesh (assume L/5 is an integer). We will also use a time step h > and consider the 
time points {ih, i = 0, 1, . . . , N}, where N — T/h. Given a mesh As, for each integer < i < N define the 



column vector F' e 



as 



/ ^0 \ 

Fl 

Fi 



V Fh J 



(54) 



with Fj = F{j6,ih). In order to measure "closeness" on this mesh, we will use the following discrete 
norm defined on As'. 

M 



U 1 



51 '^'"J I' 



u e 



(55) 



Recall that for a M x M matrix A = {aij)f^j^Q, the norm of A is given by its maximal absolute column 
sum ([5], problem 4.4.11): 



M 



\A\\i^ max V| 

7=0,1 M ^ 



(56) 



4.1 Finite-Difference Scheme for computing F{x,t) 
We begin by rewriting equation (|12p as 

V[F] = ^Fix, t) + n{t){F{x, t) - F{x - g{t),t)) = 0, 



(57) 
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where P is a linear operator defined on a suitable function space. We must first choose a discrete approxi- 
mation of V. For the time derivative, we use the usual forward difference approximation: 



, , (58) 

where the F'!- are defined in (|54p . For the difference F{x,t) — F{x — g{t),t), we use a linear interpolation 

F{xj,U) - F{x, - giU), ti)^F;-{l- X,)F;_„^ + A,i^_fe,+i, (59) 

where the integer ki satisfies Xj^^i < Xj — g{ti) < xj^ki+i, and Xi = S^^{xj — g{ti) — xj^ki)- In terms of S 
and h, these are given by 



g{ih) 



\6ki - g{ih)\ 



(60) 



See Figure [21 



\i5 



•^j — ki ■^j — ki + l 

Xj - g{ti) 



X. 



Figure 2: Relationship between Xj, giU), ki, 5 and A^. 
Putting these approximations together, we define the forward difference operator Vg : R^''^ — > M.'^^ as 

i-Pgm), = + n{ih){F; - (1 - a,)f;_,., + AO^;„,,+i) (61) 

where ki.Xi are defined in In view of ([57)1 . we require the sequence F* satisfy 

Vgm^O, i = l,2,.... (62) 

From ([HI); (|62p can we written as 

Fj^'^Fj- hn{ili)F; + h{l - X^F;_^^ + hKF;_^^^,, (63) 
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Wc can write this in matrix form. Since N is the size of time (i) mesh and M is the size of the space (j) 
mesh, we can express (|63|) as 

Y'+i = A'F\ i = l,2,...,N, (64) 
where A' is the M x M matrisQ defined as 



A' = (1 - n{7h)h)lM.M + ( Of^-Dxif-^'+i) 0(^-i)x(.-i) ^ ^ 

V ^ L>(A/-fc+l)x(A:-l) / 



(65) 



where 1m xm denotes the M x AI identity matrix, and Omxn is a M x N zero matrix, and the {M — kj + 
1) X (Af ~ ki + 1) matrix is 

{hXin{ih), 171 = 11 
h{l~X,)n{ih) 771 = n + 1 . (66) 
else 

With this, F* is easily calculated as 

F' = A'A"'K..A^F°. (67) 
where the initial condition F° is determined by the CDF of Xq and is known exactly. 



4.2 Error Analysis 

In this section, we study the error associated with solving differential-difference equation (|12p with the finite- 
difference method described in the previous section. We'll require the following assumptions on the kernel g 
and control measure n: 

• g is a monotone increasing function with Holder continuous inverse 

• g and 7i are Lipschitz functions with (68) 
\g{t) — g{s)\ < Lg\s — t\, \7i{s) — 7i{t)\ < i„|s — t\ for constants Lg, L„ > 0. 

We will see that our method gives accuracy of order o{6'^ + h^) in the discrete norm, where 7 is the 
Holder exponent oi g~^. As above, S is the size of the spacial mesh, and h is the size of the temporal mesh. 
We first consider the error associated with the approximations used in our discretization: 

dF F{xj,te + h)~F{Xj,ti) 

aF^"^"''^) " h 

F{xj - g{ti),ti) w {I - Xe)F{xj^ki,ti) + XiF{xj-ki+i,te). 

Consider the meshes = xg < xi < X2 < . . . xm and = to < ti < ••• < t^^. For integers < j < M and 
< £ < iV, define the absolute differences 

dF 

Sjj{h) = F{xj,ti + h)- F{xj,te)-h—{xj,U) (69) 
Rj,e{S) = {{l-Xi)F{xj-u,,tt) + XiF{xj^k,+i,U))-F{xj-g{U),tt). (70) 

In order to ensure convergence in our finite-difference scheme, these two errors must approach "fast 
enough" as S^h 0. Bounds for Sj^e. and Rj^i are given in Lemmas 16.11 and 16.21 in section [6] In short, these 



^Whcn implementing this method, A' should be treated as a sparse matrix to avoid running out of memory 
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Lemmas imply the following: 



\SjAh)\ 



(o{h^+^) if.T, -5(t,+i) >0 

< o{h) if Xj - g{ti,) >Q>Xj- g{ti>+i) 

[o(/i2) iixj-g{U)<Q 



For £ > 0, let F^ G M be the vector of exact values of F at time ti, i.e. 




/ F{xo,ti) \ 
F{xi,U) 



(71) 



V F{xM,te) J 



We will now state a theorem which shows that our method converges to the exact solution to ([T2l) in the 
discrete norm with a convergence rate faster than a constant multiple of [li^ + i5'*'), i.e. that ||F^ — Fg||i = 
o(/i'^ + 5'^) and ^,/i->0. 

As the proof of the following theorem will show, we will require the following stability criterion on h to 
guarantee convergence: 



where = sup n{s) and h is the time step used in the finite-difference scheme. This can always be 

0<s<T 

met easily if n is bounded. This condition makes sense since the differential equation (fT2|) holds because 
the probability of 2 or more Poisson arrivals in an infinitesimal time interval is negligible. Hence, when 
approximating (fT2|) with a finite-difference method, we must ensure that the time step we consider is small 
enough relative to the control measure to make this probability small. 

Theorem 4.1 Letg beapositive, continuous, and strictly increasing kernel with inverse g^^ G C'^''^{g{0),g{T)) 
for some < 7 < 1. Suppose the control measure n{s) is a bounded function and set — supo<f<7i 
Let Fg be defined as in ^71^ and letF^ i = 1,2, . . . , N be the approximations given by the forward difference 
scheme. Then, if hntp < 1, the forward difference scheme is convergent and satisfies 



Proof. At each step in the iteration, write A'FJ. ~ F^"*"^ + e^, where e.; G M" is the error introduced at step 
i. Then, since the initial condition F*^ is known exactly, we have 



Stability Criterion: hntj^ < 1. 



||F 



Ff 111 =o(/i^+J^) 



(72) 



pi 

F2 



F 



+eN-l+A^-hN-2 + ■■■ + A- 



(73) 
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After subtracting F"^ from both sides and taking norms, (j73p implies 

n 

Af-1 / Af-1 

< E n ii^'iii I ii^^iii 



|F^-Ff||i 



i=0 



(74) 



«=0 \'i=«+l 



where ||e£||i and ||^*||]^ are defined respectively in ([55]) and (|56p . From the definition of in ((65)l and the 
stability condition n^ft, < 1. the norm H^'Hi is given by the maximum absolute sum of the columns of A\ 
which in our case is 

Pii = |1 - nit,)h\ + \hn{U)X,\ + \hn{U)il - A,)| = 1. (75) 

Thus, dZll) becomes 



(76) 



What remains is to bound 



N-l 

\F^-Fl\\,<J2\ 

M 

\\e,\U = \\Fi+'-A'F%^Y.^\{ee),\ 

Using the differential equation (|12l) and the definitions of the errors Sj,e and Rj^i in and ([70]) . we have 

(q), ^ F{x„U+,) ~ {A'Fi), 

= F{xj,te+i) - {F{xj,ti)(l - hn{U)) + hn{te)[{l - Xe)F{xj-ke,te) + \eF{xj-k,+i,ti)]} 
= F{xj,ti + h)- F{xj,te) - h[n{te)((l - Xe)F{xj-ke,te) + XeF{xj-k,+i,ti) - F{xj,ti))] 
= F{xj,U + h) - F{xj,U) - h[n{U)(F{xj - g{U),tt) - F{xj,U) + Rj^i[5))] 

dF 



= F{xj,U + h)- F{xj,U) - h[—{xj,tt) + n{ti)Rj^i{S)] 
= Sj,i{h) - n{ti)Rj,i{5)h. 

Now, for I fixed, consider the following partition of the mesh in ([5^ : 



(77) 





= {j 


xj - g{ti>) >{)> xj- g{ti,+i)} 




= {j 


xj ~ g{te) < 0} 


4^^ 


= {.? 


X] - g{te+i) > 0} 



Notice that /-[^^ contains at most one element. Using (|77p combined with Lemmas 16.11 and 
MS = L, we have 



and that 



M 



M 



kHii = E'5K^^)^-i ^ nite)hY,s\Ra^)\+ J2 ^\s,Ah)\ + E + E ^\s,Ah)\ 

■'=0 jG/^' .e4^' ,e/<^' 

< O(W^) + • 0(5/i) + • 0(J/i2) + . 0(<5/ii+7) 

= 0(W) + 0((5/i) + 0(/i2) + 0(/ii+-^) 

= 0{h6'') + 0ih^+''). (78) 
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Finally, since Nh = T, ([751) and ([TS]) imply 



N-l 

< 0{h''+S'^), (79) 



which finishes the proof. 



4.3 Extension to more general kernels 



In order to guarantee the error bound (|72p . we must assume that g is a, strictly increasing, continuous 
function whose inverse is Holder continuous. This assumption is somewhat restrictive, as we would like to 
apply our method to a wider set of kernels. In this section, we indicate how to transform a problem with a 
general kernel into problem which satisfies the assumptions of Theorem 14.11 

Recall that if X, Y are independent non- negative random variables with respective CDFs Fx , Fy , then 
the CDF of the sum X + Y is given by the convolution 



Fx+Y (u) = Fy{u- x)Fx (dx) 
Jo 



(80) 



If Fx and Fy are both defined on a mesh = {xj}jLi, then the convolution ([80)1 can be approximated by 
the discrete convolution 

M 

Fx+y{u) « ^^^y(w - x,)[Fx{x,+,) - Fx{x,)]. (81) 

Decreasing kernels 

If g is a strictly decreasing kernel with Holder continuous inverse, then a simple variable transform re-expresses 
the integral into one with a strictly increasing kernel: 

H9)= r gis)Nids)^ rfis)N{ds) (82) 
JO Jo 

where /(s) = g{T — s) and N{ds) is a Poisson random measure with control n{s) — n{T — s). 
Flat kernels 

If 5 = A is a constant function, then there is no need to approximate the CDF since its distribution is known 
exactly: 

I{g) = I gN[ds) - APois( / n[s)ds). (83) 
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Negative kernels 

If g is a negative strietly monotone increasing/decreasing kernel with Holder continuous inverse, then the 
negative of the integral fits the proper assumptions: 

I{9) = g{s)N{ds) ^ - l^j\-g{s))N{ds)^ . (84) 

Thus, ones computes the CDF of —I{g) with the forward difference method, and then uses the relationship 

Fi(g){u) = I- F_ng){-u). 

"Piecewise" kernels 

Combining the above three ideas, the convolution formula (|8ip and the independent increment property 
of the random measure N ^ we can approximate integrals with piecewise increasing/flat/decreasing, posi- 
tive/negative kernels, whose inverse on each non-flat piece is Holder continuous. Indeed, if g is increas- 
ing/flat/decreasing and positive/negative between each of the points = to < < ■ ■ • < tn-i < tn = T, 
then we have 



T 



i+i 



I{g) = / g{s)N{s) ^ ^ / g{s)N{ds) (85) 

This sum on the right-hand side can be computed inductively by using the integral up to ti-i as an initial 
condition, and iterating until the CDF is computed up to ti. That is, set Jq = 0, and use our method to 
compute 

/,=/,-!+/ g{s)N{ds), ^ = l,2,...,n. (86) 
•Iti-i 

Corollarv 13 . 2 1 and Thcorcm l4. II implv that the rate of convergence in this method in the norm will depend 
on, at worse, the minimum of the Holder exponents of g~^ over each of the intervals [t^, i^+i). 



5 Examples 

We shall demonstrate our method on some examples. We start with 

/ sN{ds), (87) 

where iV is a Poisson random measure with control Lebesgue: i.e. n(s) — 1. This example is simple enough 
to compute exactly, so it serves as an useful test case. The second example is where the kernel g is a parabola 

\l-{l-sf)N{ds), (88) 

and the control measure is again Lebesgue. In this example, the kernel g has an inverse which lies in the 
Holder space C°'^/^(0, 1), and thus we expect the CDF to share this property. For purposes of comparison, 
we also compute this CDF by approximating the integral ([5]). 

Finally, we consider an integrand which is both positive and negative: 

/ sin(27rs)iV(ds) (89) 
Jo 
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with control measure Lebesgue. This will demonstrate some of the ideas in Section [ 

Recall that in each example, we will see a jump at a; = of size exp(— n(s)ds) in the CDF which is 
the probability of no Poisson arrivals. 



Example 1: 



sN{ds), n{s) = 1 



In this example, the kernel is the identity function ^(.s) ~ s, hence the integral is given by the sum of arrival 
times in (0,1). The number of arrivals follows a Poisson distribution with rate n{s)ds — 1 and in the 
event of fc > arrivals, the times of these arrivals are given by k i.i.d. random variables Ui, . . . ,Uk which 
are uniform on (0, 1). Thus, the CDF F of Jg sN{ds) can be expressed by conditioning on the number of 
Poisson arrivals in the interval (0, 1): 



F{x) = ^ Prob 



A:=0 



i=l 



k arrivals 



Prob[fc arrivals] = 



fe=0 



k\ ' 



(90) 



where Pk{x) is the CDF of a sum of k i.i.d. C/(0, 1) random variables. Each Pk can written as a piecewise 
polynomial, and can be computed recursively as 



^^0(2;) 
Pk{x) 

The first few Pt are 




/ Pk^i{x — u)du k > 1. 
Jo 



Piix) 



0, x<0 
X, < a: < 1 

1, x>l 



P2{X) 



'0, 



2 ' 

(-2+4a: 



X < 
< X < 1 

, Kx <2 



Psix) 



.1, 



2 

x>2 



0, X < 

^, < a; < 1 

3-9a;+9a:^-2x^ 



-21+27a:-9a;^+a^ 
6 

1, X > 3 



1 < X < 2 
2 < X < 3 



Notice that it suffices to take the first 11 terms in the sum (|90p . since X^i^ii e^^/n! < 10^ , which is much 
more than the precision we seek. 

We generated F{x) on the interval x G [0, 3) with our method with S = h = 10^^. This is plotted in 
Figure [3] along with a numerical approximation of the densitjU. Notice the "kink" at x = 1 in the CDF and 
the corresponding discontinuity in the density. This is consistent with Theorem 13.11 since here the inverse 
of the kernel is simply g~^{s) = s, which lies in the space C'^'^(0, 1), so here we take 7 = 1. Thus, we expect 
that the F is at worst Lipschitz continuous. This is indeed the case, since P is a linear combination of the 
Pji's, and Pi is clearly in C°'^(0, 00) and no smaller Holder space. 

Using the first 11 terms in the sum (|90|) . we also computed this CDF on the interval [0, 3) and looked at 
the relative error between this and the result of our method. The relative error is plotted in Figure S) We 
have excellent agreement over all x with our method. 



This is obtained using the central difference F'(x) Ri (25i)-i(F(x + 5i) 
for the error made in computing the CDF with step size 5. 



F{x — Si)). We must choose <5i > <5 to account 
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Figure 3: CDF and density for the stochastic integral ([87|) . The vertical line with a dot in the plot of 
the density represents a discrete atom in the distribution whose mass is the height of the vertical line. To 
compute this, we used the forward difference method with h = S = 5 x 10^**. On a 1.5 Ghz computer, this 
computation took less than 5 seconds using MATLAB. Notice the "kink" in the CDF and discontinuity in 
the density at a; = 1. This is not surprising, since here 7 = 1 so expect the CDF to be at worst Lipschitz 
continuous. 

Example 2: [ [I - {1 - sf]N{ds), n{s) = I 
Jo 

Here we consider the more complicated example where g{s) — 1 — (1 — s)^ for < s < 2, which is a downward 
facing parabola with maximum value 1 and roots at and 2. This function is piecewise monotone on the 
intervals (0, 1), and (1, 2), so the finite difference scheme can be applied first to compute F{x, t) {01 < t < 1, 
and then subsequently for 1 < i < 2 by using F(x^ 1) as an initial condition. The density can be obtained 
with numerical differentiation. This is shown in Figure \5\ 

Notice the sharp corner in the CDF at .t = 1 and the corresponding singularity in the density. This 
illustrates the results of Theorem 13.11 and Corollary 13.31 that is, the CDF lies in the smallest Holder space 
C^''^{0, 00) which contains the inverse of g restricted to (0, 1) and the inverse of g restricted to (1, 2). In this 
case we have 

g^^{x) = 1 — \/l — x, for g restricted to (0, 1) 
g~^{x) = 1 + \/l — X, for g restricted to (1,2) 

In both of these cases, g^^ £ C^'^^^{0, 1), and thus we expect the CDF F{x,t) e C°-^/^{0,oo) as well. 

For comparison purposes, we also computed this CDF using the integration method of Abate and Whitt 
discussed briefly in the introduction. Since this integral is only applicable for continuous F, we first rewrite 
the characteristic function (j)g by "removing" the case where there are arrivals, and considering this case 
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10 ■ 



Relative error 




Figure 4: The relative error between the CDF of ([57]) computed using our method, and the true value of 
the CDF computed with (PO)) . We see that our method gives accuracy of ^ .1% over all x considered. 



separately, 



Eexp ( ie J g{s)N{ds) 



Prob[A^ = 0]e° + E 
r2 



exp [ I g{s)N{ds) 



N > 1 



exp [ / (e 



iOais) 



- l)ds - e 



Then, using ([5]) the CDF F is given by 

F{x) = e-^ 



Re((Pg)(M) du. 



(91) 



To approximate the integral in (|9ip . we truncate the infinite limit at T/ > 0, choose a mesh size > and 
use a trapezoid approximation: 



Fix) « + I ('Re(0,)(O) +Re(^,)(T,)^^^^ 



(92) 



Where X = Tj/rj is the size of the mesh. Notice that computing Rc{(j)g){jr]) for every j involves also 
computing the integral /q^ (e'^^'-"^ — l)ds. To do so, we used the MATLAB function quad. The integrand is 
a highly oscillating function for large 9, causing numerical approximation methods to converge slowly. This 
makes computing the sum in (|92[) extremely time-consuming for large K. 

In Figure [51 we plotted the CDF around the kink at x = 1 using these two methods. We used the 
finite-difference method with S = h = 5 x 10^^ and S = h = 1 x 10^^, and used the integration method 
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with T/ = 50, r/ = .1 and T/ = 100 and ry = .01. Notice that there is a small difference between the 
curves generated by the finite-difference method, which implies adequate convergence has been met. On the 
other hand, the results of the integration have "smoothed" out the kink, implying a larger T/ is required. 
Moreover, the computation time for these methods varied drastically, to generate the curves in Figure [SI the 
finite-difference method required less than 2 seconds, while integration requires about 8 minutes. 



CDF 



Density 





Figure 5: CDF and density the Poisson integral of g. In this case, the CDF is Holder continuous for x > 
with Holder exponent 7 = 1/2. As a result, we see a sharp "corner" in the CDF and a singularity in the 
density at a; = 1. To generate the CDF, we used 6 = h = 10^*. 



Example 3: f sm{2TT s)N{ds), n{s) = 1 
Jo 



We now consider the Poisson integral of one period of a sine wave which is a case where the kernel is both 
positive and negative. To compute the CDF F, we split up the integral as 

sm{2TT s) N (ds) ^ [\in{2TTs)N (ds) + [ sm{2Tr s) N (ds) =: h + I2 (93) 
Jo 

and note that the summands on the right are independent and satisfy Ii = —l2- Thus, if Fi{x) is the CDF of 
/i, then -^2(2;) = 1 — Fi{—x) gives the CDF of l2- Taking the convolution of Fi and F2 we use the symmetry 
of Ii and I2 and see 



F{u) = I F^{u - x)dF2{x) 

-00 

poo 

■■^/^Fi{u)+ / Fi{u + x)dFi{x) 



e 







M 

e-^'^F^{u) + Y. ^1 + ) (^1 (^J- )-Fi{xj.i)). (94) 
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Figure 6: A comparison between the CDF of Example 2 computed using the integration method (|92|) and 
the finite-difference scheme of Sectional zoomed in on the kink in the CDF at a; = 1. The integration 
method takes over 200 times longer to implement, and it seems that even larger Tj and smaller j] is required 
to obtain sufficient convergence. Also, important features of F are missed with integration, such as the kink 
at a; = 1 which is "smoothed" here. The finite-difference method, on the other hand, converges quickly in 5 
and h, is much faster, and captures the kink in F 

The first term appears here since F2 has a discontinuity at a; = with size 

F2(0+) -F2(0-) = e-^o^'^id',) ^ 

This gives an approximation to the CDF of the sum in ((93|) . 

The CDF and density of sin(27rs)iV(ds) is plotted in Figure[7l Similar to Example 2, we see at a; = ±1 
two sharp corners in the CDF. To understand this, notice that on the interval (0, 1/2) the inverses of the 
function g(s) = sin(27rs) are given by 

Sin f 

g^^{x) = — , for 5 restricted to (0, 1/4), 
sin (1/2 — x) 

g'^ix) = Y ^ for g restricted to (1/4, 1/2). 

In both of these cases, g^^ G C°'^/^(0, 1) since the derivative at x = 1 doesn't exist and 

sin^"'"(l) — sin^^(a::) ~ \/l — x, asx^l^. (95) 

Thus, Theorem [O] and Corollarv l3 . 31 imply that Fi G C°'^/^(0, 00), and similar to Example 2, this is evident 
by a kink in in Fi at a: = 1. Since F2{x) = 1 — Fi{—x), F2 G C"'^/^(— 00, 0), and we expect a kink at a; = — 1. 
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Figure 7: CDF and density the Poisson integral of the sine wave g{s) = sin(27rs) over the interval < s < 1. 



Because Fi and F2 both have a discontinuity at a; = 0, the convolution in ()94|) can be written in two ways: 



Fix) 



Fi{x-u)dF2{u)^e-^^^Fi{x)+ / Fi{x - u)dF2{u) 

> J —00 

F2{x-u)dFi{u)^e-^'^F2{x)+ / F2{x - u)dFi{u). 



(96a) 
(96b) 



Since Fi has a kink at x = 1 and F2 has a kink at x = —1, ()96ap and (|96bp explain why you see kinks at 
a; = ±1 in the convolution. 



6 Approximation Lemmas 

This section contains the bounds on the approximation errors Sj^i{h) and Rjj{5) defined in Scction|4?2l Recall 
that Sjj{h) involves an approximation of the time derivative and Rj^i{5) involves an approximation of 
F(x — g(t), t). The first lemma below bounds Sj,i{h) in terms of h, h^^'' and . The second lemma bounds 
Rj^i{6) in terms of 5^+'^. Here 7 is the Holder exponent of g^^. 

Lemma 6.1 Let Sj^i{h) he as in ^69jl. and assume g and n satisfy the assumptions in Ii68\} . Then 



n^{2h + n^h ) Xj - g{ti) > > - g{te.+i) 

Z 1 + 7 

^ + 2n1^)h'^ Xj - g{tt) < 

where Cn,T is defined in Theorem \S.l\ 



(97) 
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Proof. Let h > 0. As in the proof of Theorem 12.11 we have 

F{xj,te + h)- F{xj,U) = E[l(_„c..](^(t£ + M)-l(-oo,.](^(i£))] 

OO 

= ^E[l(_oo,.](^(t« + h)) - l(_^,,](X(t,))|iV = 7i]P„ 



n=0 



+ IE[l(-oo,.] + h)} - l(_oo,.] = n]Pn 

Now, using dni, dMl), and the fact that < F(xj,tt) < 1, we have 

OF 

\S,Ah)\ = \F{xjM + h)^F{x,M)-h—\ 

= \F{xj,ti + h)- F{xj,U) - h{F{xj - g{te)M) - F{xj,ti))n{U)\ 

ti + h \ 



(98) 



< 



(F{xj ~ g{s)M) - F{xj,h))n{s)ds - h{F{xj - g{U),tt) - F{xjM))n{ti) 



ti+h 



{F{x, - g{s),U) - F{x,,U))n{s)ds) (g" _ i) 



n=2 



(99) 



To obtain ()97p . we consider the three cases separately. We will use the following elementary inequalities 
involving e^^ for x G K: 



l-e-'^ < X 
1 - e^"" - xe^"" < x^ 



(100a) 
(100b) 



Case 1. Xj — g{te) > > Xj — g{t(+i) 

This case contains the "troublesome" point, since ^ C{F) and Theorem 12.11 does not necessarily apply. 
However, in view of (jl8p in Remark 2 after the proof of Theorem 12. 1[ we have in any case 



dF 

I r 

dtt 



< hrin 



(101) 



And, from 



\F{xj,U + h)-F{xj,ti)\ < n^^h + J^Pn- 



n=2 



Using (jlOObp . we have 

OO 



t+h 



< 



t+h 



n{s)ds < (/in^)^. 



(102) 
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Thus, for this case, 



dF 



< 2n*rph 



(103) 



\SjAh)\ < \F{xj,U + h)- F{xj,U)\ + 

which is consistent with (j97|) . 
Case 2: xj — gitg+i) > 0. 

Since g is monotone increasing, Xj — g{ti) > Xj — g{t£-i-i) > 0. From Theorem 13.11 and the Lipschitz 
assumptions on n and g imply the first term in absolute values in (|99p can be rewritten and bounded as 



rti + h 



{F{xj - g{s), ti) - F{xj - g(ti),ti))n{s)ds 



tt+h 



{n{U) - n{s)){F{xj - - F{xj,ti))ds 



tt+h 



\F{xj - g{s),ti) - F{xj - g{ti),te))\ds + / \n{s) - n{ti)\ds. 



tt+h 



tt+h 



< n*j.Cn,T[g ^]co.Ug{o),g(T}} \g{s) - g{t)\'^ds+ \n{s) - n{ti)\d. 



ti 

tt+h 



tt+h 



ntt+h 

\s — tepds + Ln / \s — te\ds 
Jt, 



= n*rpCn,T[g ^]c0'-'{a{0),g[T))LlY^+ ^^^^ 
From (|100ap and since < |^| < 1, the second term in ((99)) is bounded by 

ftt+h 



tt 



< hntp 



tt+h 



n{s)ds 



[Fix, - g{s),U) - F{x,,U))n{s)ds^ (g- " _ i) 

Finally, the third term in is bounded by (fT02| . Putting (fT04|) . (fTOS]) and (fTOS]) together, gives 



\Sj,i{h)\<{— + 2)nT, h + — — 

z 1 + 7 



(g(o),g(t))-^g ^i+7 



(104) 



(105) 



(106) 



which agrees with (^7]) for this case. 
Case 3. Xj — g{ti) < 

In this case, many of the terms in (|M)) drop out since F{x, t) = for x < and g is assumed to be 
non- negative. Breaking up (|99| in a similar way to case 2, we obtain 



\SjAh)\ < 



tt+h 



\n{s) -n{te)\F{x.j,tt)ds + 2n*j.'^h'^ 



< (!!z:il!i + 2n5,)/i2 



This confirms ([57)1 and finishes the proof. 
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Lemma 6.2 Let Rj^^{5) he as in \7(]^ , and assume g satisfies the assumptions in i [6'<^|) . Then 

\R,,im < C„,t[.9"']co,.(3(o).s(t))<5^ (107) 
where Cn.T is defined in Theorem \3.1\ 

Proof. By the definitions of and kg, and since F and g are non-decreasing, we have 

F{xj^ki,ti) < F{xj ~ g{tt),ti) < F{xj-kt+i,ti) 

F{xj-ke,te) < {I - Xe)F{xj-ki,ti) + XeF{xj-ke+i,ti) < F{xj-ke+i,t£). 
These inequalities and Theorem 13. II now imply 



\RjA^)\ ^ \F{xj_ki,+i,ti) - F{xj_ki,ti)\ 

< Cn,T[g~^]co.-'{g{0),g{T))\Xj-kt + l " Xj- 
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= CnMg '^]co.-'{g{0),g{T))S'^ ■ (108) 

This finishes the proof 



7 Guide to Software 

Software written in MATLAB for computing the CDF of a Poisson integral of the form ([1]) is freely available 
from the authors. This section will go over installation and use. 

Once the file Poisson_Integral_GUI . zip is downloaded and decompressed, move the folder "Poisson_Integral_GUI" 
to your desired directory. Launch MATLAB, and enter path (path, 'your4)ath/Poisson_Integral_GUI ' ) , 
where your_path denotes the path to the folder Poisson_Iiitegral_GUI. You should now be able to use the 
code included with the package. 

To begin, type Poisson_Integral_GUI and press return. This will open a window similar to that shown 
in figure m From here you may enter the parameters into the "Inputs" panel on the top left. These are: 

5 Step size in spatial mesh 

h Step size in temporal mesh 

T Upper limit of integration in ([1]) 

Xmax Largest value for which CDF is computed 

g{s) Kernel function 

n{s) Density of the control measure 

The kernel function and control measure should be entered as MATLAB expressions in terms of s. For 
example, the kernel g{s) = sin^(27rs) and Lebesgue control measure n(s) = 1 are entered as sin (2 * pi * 
s)A2 and 1, respectively. Only non- negative kernels and control measures will be accepted. 

Once all the parameters are entered, clicking on Compute F will generate the CDF of the integral 
Jo 9is)N{ds) in the upper right plot, as well as a "quick and dirty" approximation of the continuous part 
of the density in the lower right plol|j. The density is approximated by fitting a smoothing spline to the 



^This requires the Spline Toolbox. 
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Figure 8: A snapshot of the GUI. 



CDF. then computing the derivative of this spUne. The shder under the density plot aUows you to adjust 
the smoothing parameter used in computation of the density. 

Once a CDF is computed, you may compute specific values of the CDF using the "Outputs" panel on the 
left. Simply enter in a desired value in the x slot, and press return. Alternatively, you may select an a;-value 
using the slide bar at the bottom of the panel. You can also enter a value in the F{x) slot to compute the 
inverse of the CDF. 

The lower left panel allows you to save the CDF data to your MATLAB workspace. To do this, enter a 
variable name, and click save. This will create a variable on your workspace with your chosen name. This 
variable is an M by 2 array whose first column contains the x values (0, 6,26, ... , M6y and second column 
contains the corresponding CDF values, i.e. {F{0),F{6), . . . , F{M6))' (here, M denotes the size of the spatial 
mesh). 

Finally, if the user wants to avoid using the interface, the code poissonCDF will work in the command 
line. This function takes in the six parameters listed above and outputs a vector containing the values of 
the CDF on the mesh points generated. The functions g and n must be entered as function handles. For 
example, \i 6 = 0.01, h — 0.01, T = 3, Xmax = 4, g{s) = and n(s) = 1, the call to this function would look 
like 

F = poissonCDFC .01, .01,4,3,@(s) s A 2,@(s) 1 ); 

The "0" symbol is necessary because the input must be a so-called "function handle" . The statement above 
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will return the vector F which contains the values of the CDF of the integral s'^N{ds) with Lebesgue 
control measure for the values x — .01, .02, . . . , 3.99, 4. For more on this function, read the comments in the 
code. 
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